//////////////////////////////////////////////////////////////////////
// sparse_binding.h
//
// SHEN Fangyang
// me@shenfy.com
//
// Copyright (C) SHEN Fangyang, 2011, All rights reserved.
//////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////
//  Brief Description
//
//	SIMPLE LINEAR ALGEBRA PACKAGE
//
//	Binding CSparseMatrix and ublas vectors to Intel(R) MKL
//	including sparse BLAS, LAPACK
//
//////////////////////////////////////////////////////////////////////
#pragma once

#include "sparse_detail.h"

namespace grflib
{
	namespace LA
	{
		namespace bindings
		{
			//Sparse matrix-dense vector multiplication
			//y <- alpha * (A^T | A) * x + beta * y;
			template <typename matrix_type, typename vector_type_x, typename vector_type_y>
			inline void csrmv(const char TRANS, //'T' or 't' to transpose, 'N' or 'n' if not
				const typename matrix_type::value_type &alpha,
				const matrix_type &a, //must be a CSparseMatrix
				const vector_type_x &x,
				const typename matrix_type::value_type &beta,
				vector_type_y &y)
			{
				using namespace boost::numeric::bindings;

				const int m = static_cast<int>(a.m_meta.nRow);
				const int n = static_cast<int>(a.m_meta.nCol);

				//const int stride_x = traits::vector_stride(x);
				//const int stride_y = traits::vector_stride(y);

				const typename matrix_type::value_type *a_ptr_value = a.m_values;
				const int *a_ptr_indices = (const int *)(a.m_indices);
				const int *a_ptr_B = (const int *)(a.m_ptrB);
				const int *a_ptr_E = (const int *)(a.m_ptrE);
				
				const typename matrix_type::value_type *x_ptr = 
					boost::numeric::bindings::traits::vector_storage(x);
				typename matrix_type::value_type *y_ptr = 
					boost::numeric::bindings::traits::vector_storage(y);

				detail::csrmv('N', 
					m, n, 
					alpha, 
					"G__C", //general matrix, zero-based
					a_ptr_value, 
					a_ptr_indices, 
					a_ptr_B, 
					a_ptr_E,
					x_ptr,
					beta,
					y_ptr);	
			}

		} //namespace bindings
	} //namespace LA
} //namespace grflib